#include <stdafx.h>
#include <iostream>
#include <fstream>
#include <stddef.h>
#include <stdexcept>
#include "blinalg.hpp"
#include "usertypes.hpp"
#include "string.h"
#include "gauss.hpp"
#include "superelement.hpp"
#include "functions.hpp"
#include "Intergration.hpp"

double Intergration(double (*f)(double, FElement) , const FElement elem, int nr_gauss)
{
	double x, d = 0;
	init_gauss();
	double a = elem.getVertexCoord(0);
	double b = elem.getVertexCoord(1);
	//
	for(int i = 1; i <= nr_gauss; i ++)
	{
		x = a+(b-a)*xi[i][nr_gauss];
		d += wi[i][nr_gauss]*f(x, elem);
	}
	d *= elem.getArea();
	return d;
}
double Intergration(double (*f)(double, FElement,EVector<double>,vector<FElement>) , const FElement elem, int nr_gauss, EVector<double> R, vector<FElement> m_elems)
{
	double x, d = 0;
	init_gauss();
	double a = elem.getVertexCoord(0);
	double b = elem.getVertexCoord(1);
	//
	for(int i = 1; i <= nr_gauss; i ++)
	{
		x = a+(b-a)*xi[i][nr_gauss];
		d += wi[i][nr_gauss]*f(x, elem, R, m_elems);
	}
	d *= elem.getArea();
	return d;
}

double Intergration(double (*f)(double) , int nr_gauss)
{
	double x, d = 0;
	init_gauss();
	double a = 0;
	double b = 1;
	//
	for(int i = 1; i <= nr_gauss; i ++)
	{
		x = a+(b-a)*xi[i][nr_gauss];
		d += wi[i][nr_gauss]*f(x);
	}
	return d;
}
double Bii_1(double e, FElement ip_elem)
{
	double a = ip_elem.getVertexCoord(0), b = ip_elem.getVertexCoord(1);
	double h = b - a;
	//double x = h*e + a;
	return Function_p<double>(e)/(h*h) + Function_q<double>(e)*((e-a)/h)*((e-a)/h);
}
double Bii_2(double e, FElement ip_elem)
{
	double a = ip_elem.getVertexCoord(0), b = ip_elem.getVertexCoord(1);
	double h = b - a;
	//double x = h*e + a;
	return Function_p<double>(e)/(h*h) + Function_q<double>(e)*((b-e)/h)*((b-e)/h);
}
double Bii1(double e, FElement ip_elem)
{
	double a = ip_elem.getVertexCoord(0), b = ip_elem.getVertexCoord(1);
	double h = b - a;
	//double x = h*e + a;
	return (-1)*Function_p<double>(e)/(h*h) + Function_q<double>(e)*((b-e)/h)*((e-a)/h);
}
double F_1(double e, FElement ip_elem)
{
	double a = ip_elem.getVertexCoord(0), b = ip_elem.getVertexCoord(1);
	double h = b - a;
	double x = h*e + a;
	return Function_f<double>(e)*(e - a)/h;
}
double F_2(double e, FElement ip_elem)
{
	double a = ip_elem.getVertexCoord(0), b = ip_elem.getVertexCoord(1);
	double h = b - a;
	double x = h*e + a;
	return Function_f<double>(e)*(b-e)/h;
}
double F_error(double e, FElement ip_elem, EVector<double> R, vector<FElement> m_elems)
{
	//double a = ip_elem.getVertexCoord(0), b = ip_elem.getVertexCoord(1);
	//double h = b - a;
	//double x = h*e + a;
	return (ExactSolution<double>(e) - app_solution(R, m_elems, e))*(ExactSolution<double>(e) - app_solution(R, m_elems, e));
}
double app_solution(EVector<double> ip_c, vector<FElement> ip_elems, double ip_x)
{
	int v_n = ip_c.size(), i0;
	double v_tmp = 0;
	i0 = findIndex(ip_elems, ip_x, 0, v_n-1);
	if (i0<v_n-1)
		v_tmp += ip_c[i0] * (ip_x - ip_elems[i0].getVertexCoord(0))/ip_elems[i0].getArea();
	if (i0>0)
		v_tmp += ip_c[i0-1] * (ip_elems[i0].getVertexCoord(1) - ip_x)/ip_elems[i0].getArea();
	/*for (int i = 0; i < v_n; i++)
	{
		if ((ip_x > ip_elems[i].getVertexCoord(0)) && (ip_x < ip_elems[i].getVertexCoord(1)))
			v_tmp += ip_c[i] * (ip_x - ip_elems[i].getVertexCoord(0))/ip_elems[i].getArea();
		else if ((ip_x > ip_elems[i+1].getVertexCoord(0)) && (ip_x < ip_elems[i+1].getVertexCoord(1)))
			v_tmp += ip_c[i] * (ip_elems[i+1].getVertexCoord(1) - ip_x)/ip_elems[i+1].getArea();
	}*/
	return v_tmp;
}

int findIndex(vector<FElement> ip_elems, double x, int d, int c)
{
	if (d>c) return d;
	int g = (d+c)/2;
	if ((ip_elems[g].getVertexCoord(0) <= x) && (ip_elems[g].getVertexCoord(1) > x))
		return g;	
	if (ip_elems[g].getVertexCoord(1)==x)
		return g+1;
	if (ip_elems[g].getVertexCoord(0)>x)
		return findIndex(ip_elems, x, d, g - 1);
	if (ip_elems[g].getVertexCoord(1)<x)
		return findIndex(ip_elems, x, g+1, c);
	return d;
}